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FOREWORD 

This document presents the results of work performed by 
the Fluid Mechanics Section of the Aeromechanics Department 
of Lockheed's Huntsville Research & Engineering Center, This 
report is Volume III of a four -part final report, as required to 
fulfill Contract NAS7-761. This work was sponsored by the 
Liquid Propulsion Section of Jet Propulsion Laboratories, Mr. 
Wolfgang Simon, Technical Manager. 

This document constitutes Volume III of the four -part 
final report. The other three volumes, printed separately 
are: 

Volume I — "Summary Volume — Method of Character- 
istics Nozzle and Plume Programs," LMSC-HREC 
D162220-I. 

Volume II — "User's Manual — Method of Cha»acter- 
istics M.ot Program, ' LMSC-HREC D162220-II. 

Volume IV — "User's Manual — Variable O/F Method- 
of-Characteristics Program for Nozzle and Plume 
Analyses," LMSC-HREC D162220-IV. 


..EMSION NOTICE 

Lockheed Missiles & Space Company, Huntsville Research 
& Engineering Center Technical Report HREC 7761-3 (D162220-III), 
entitled "Solution o|f Non-Isoene- getic Supersonic Flows by Method 
of Characteristics," dated July 1971 is revised as indicated below: 


Revision A affects title page and pages ii, iii and 68. 
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SUMMARY 


This is one of a series of reports dealing with the calculation of 
supersonic flow fields by the method of characteristics. This report deals 
with the theoretical approach to the solution of these flow fields while sub- 
sequent reports will compare the theory w'ith experiment and discuss in 
detail a computer program which was derived to implement the numerical 
solution of the flow equations. This versatile program has a flexible set 
of boundary conditions enabling the calculation of nozzles, plumes and many 
other complex flow fields. A user's guide for this program is contained in 
Volume IV, "User's Manual -* Variable O/F Method-of-Characieristics Pro- 
gram for Nozzle and Plume Analyses," LMSC-HREC D 162220 -IV. 


A complete derivation of the equations of motion for re.>cting gas 
systems is presented in this report. This derivation clearly illustrates the 
underlying assumptions that were made to arrive at the more familiar sys- 
tem of equations which was finally treated. An important consequence of 
this derivation is that, for the reaction assumptions which were made, the 
thermochemistry was shown to be uncoupled from the flow solution and as 
such could be solved separately. In addition the method of characteristics 
equations are shown to be formally the same for ideal, frozen, and equili- 
brium reacting gas mixtures. 


The two dimensional and axisymmetric characteristic equations are 
cast in finite difference form. These equations' apply only in regions of the 
flow field where transport properties can be ne|glected. A shock wave solu- 

‘ t 

tion is pr<isented, and numerical techniques necessary to effect the calcu- 
lation are described. The shock wave solution, then acts as a patching line 

between r igions in which characteristic re^tions ^pply. 

t i 

1 ! ji 

A di Bcussion of the mesh construction nec^s^ary to apply the local 
point solutions in the proper fashion in also givle 

: Jill: 

iii 
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SYMBOLS AND NOTATIONS 
Description 

time 

volume 

area, chemical symbol 
mass density 
velocity 

diffusion velocity 
reaction coefficient 
mass 

progress variable 
stress tensor 
unit tensor 

viscous stress tensor 
body force vector 
hydrostatic pressure 
interi^l energy per unit mass 

heat flux vector ii ’ 

1 ! 

enthalpy per unit m^ss 
entropy per unit mass ; 
temperature . ! . 

mass fraction ) H 

r » H 


1 
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Symbol 
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X® 

W 

y 

R 

2 

a 

u, V 
X, T 

e 

n, t 
8 
€ 

M 

E 

<1> 

j.k 

V 

D/Dt 


Description 

Gibbs potential, Mach angle 

Gibbs free energy change 

equilibrium progress variable 

universal gas constant 

molecular weight 

isentropic exponent 

local value R /W 
u m 

speed of sound squared 
velocity components 
characteristic slope 
position coordinates 

inclination of flow vector with respect to x axis 

normal and streamline direction 

equation modifier, turning angle 

shock angle 

Mach number 

error function 

necessary information to describe point 
operation function 


known and partially known characteristic 


Nabla 




i: 


substantial derivative* 

■ i i r 
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Superscripts 
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R 
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Description 

vector, also denotes averaged quantity 
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N. E. 


i species, otherwise where defined 

k reaction 

summation index 

pertaining to the mixture 

equilibrium and frozen 
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Section 1 
INTRODUCTION 


The method of characteristics is known to be an accurate means of 
calculating supersonic flow fields. Since the advent of modern high speed 
digital computers, this method has recieved much attention because of its 
accuracy and, in spite of, its complexities. The equations found in this 
report, although perhaps different in minor ways, are those which have 
been known and used for half a century. It is felt though that the manner 
in which these concepts have been combined to result in a useful, flexible 
tool for rapid calculations is unique. 

A complete derivation of the equations of motion for reacting gas' 
mixtures is presented in this report. This was done so that the conse- 
quences of assumptions necessary to arrive at the simple forms ulti- 
mately treated can be examined. For the special types of reactions 
treated it is shown that the thermochemistry can be completely separated 
from the flow solution thus greatly simplifying the entire calculational 
procedure. 

In spite of the almost universal application which the method has 
recieved, very f^w descriptions of the mesh construction techniques 
necessary to apply this method have been documented. It is the intention 
of this report to present this important portion of the theory. A functional 
description of the mesh construction will be attempted in the hope that it 
will enhance the understanding of this facet of the technique. 
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Section Z 

THE FUNDAMENTAL EQUATIONS FOR FLOW 
OF REACTING GAS MIXTURES 


2.1 DERIVATION OF EQUATIONS FOR UNSTEADY NON EQUILIBRIUM 
FLOW 


2.1.1 The Continuity Equation 

In a flow system of gaseous mixtures in which chemical reactions 
take place, the principle of conservation of mass of each chemical species, 
when applied to a control volume V bounded by a control surface A, may 
be written as 

V A V 


1 ^ X I • • • • A 

th 

where p. is the partial density of the i species; ^ the ^^elocity vector of 

^ th ^ 

the center of mass of the i^^ speciPo> dA the utward differential area 

vector on the control surface A; dV the differential voliuno element; 
the rate of increase of due to either internal or external sources such as 
chemical reaction and mass additions; and n is the total number of chemi- 
cal species in the mixture. I 


th 

Th«i above equation simply states that, fpr the i species, the masi 
rate of increase inside the control volume is ^qual to the net rate of mass 
flow into the control volume plus the rate of mass produced due to chemi- 
cal react ions and mass addition. In the following, we shall, however, 
exclude riass addition from external sour<^es.; 




I 
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The first integral on the right-hand side of Equation (2-1) may be 
transformed by means of Gauss^ theorem 


i 


t) • dA 


V* 5(r.t) dV 
V 


where i? (r, t) may be any continuous vector field (functions of position 
vector r and time t) defined inside and on the boundary of the control vol- 
ume V. Also, in the Eulerian coordinate system which we use here, the 
control volume is not a function of time, and the order of differentiation 
and integration may be interchanged, i.e.. 


V 


dV 



Hence, Equation (2-1) becomes 



( 2 - 2 ) 


Since the volume V under consideration is arbitrary, the only way for 
Equation (2-2) to be valid for all V is for the integrand to vanish. We 
therefore have 


®p. 

+v [Pi(q + ^o] 


hi 

W 


(2-3) 


i = 1.2, 


n 


respect 


where wi have replaced q. by q + 5^. q is the| velocity of the center of mas* 


of the mixture at (f.t) and is the diffusion velocity of the i^^ species with 


o the mass center of the mixture^ 


I I 


3 


I 
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Equation (2-3) is known as the species continuity equation. It is 
valid for each^chennical species at each internal quantum state. \Vc shall, 
however, assume that the various internal modes of motion are fuUy ex- 
cited and are in equilibrium with each other. It is well known that this is 
approximately the case for the translational and rotational degre^?s of 
dom where the equilibrium value is attained in a few collisions. In general, 
the vibrational degree of freedom approaches the equilibx-iuni state some- 
w'hat more slowly, except at very high temperatures. As cheriical reac:ions 
usually occur at high temperatures, this approximation is of‘en justified. 


If one sums up the n equations represented by Equation (2-3) a.'d 
utilizes the relationship between the mass density p of the mixture and the 
species density* p. 

n 

P = (2-4) 

i=l 


the requirement for the diffusion velocity that 



pq - pq = 0, 


(2-5, 


and the conservation of total mass for chemical reactions. 


one readily obtains 


t 

i=l 


Sp, 

It 


= 0 




^ + V • (pq) = 


( 


(’-7) 


Eqfiation (2>7) has the usual fotm of ma^s conservation or global 
continuit|r equation. It c\n be written alternatlrefy as 

' 'nil 
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•D£ 

Dt 


+ P ^ • q 


0 


(2-7a) 


by means of Euler's derivative 


Dt ■ dt 


+ q .7 


The species continuity Equation (2-3) may be written in more con- 
venient forms. By making use of the global continuity equation and intro- 
ducing the mass fraction c.(=p./p), Equation (2-3) becomes 


Dc. 

+V-(PCiqp 



(2-8) 


The term on the right-hand side of Equation (2-8) may be put into 
terms of progress variables of the chemical reactions. 

Supposing tl. jre are r independent reactions occurring in a mixture of 
n species, then a typical reaction may be represented by the general 
form 


n 

Aj * 0 , k si, 2 , ,,,• r 
U1 


(2-9) 


where is the chemical symbol for component i and is the number of 
grams of component i produced per gram of reaction U. Since the mass 


of a re 


actingj 


system is not, it is more convenient to use the gram basis rather than the 


mole basis i^ 
is negative ii 


system is conserved while the number of moles in such a 


h ou~ discussion of chemical reac'.ion^. The coefficient V. 


ik 

i is a reactant i*.nd positive if i is a product of the reaction k. 


S 
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If component i does not participate in reactions k, then the corresponding 

V.. is sero. 
tk 

It follows, obviously, from the above description that if there is no 
mass addition from external sources, 

2 = 0. k = 1.2 r (2-10) 

U1 


Furthermore, the changes of mass of the components in the system 
th 

as a result of the k reaction are related by 


^•lk 




2k 


d. m 

, k = 1. 2 r (2-11) 

*^nk 


Now we introduce a progress variable for the k^^ reaction and 
define it as the grams of reaction occurred per gram of original reactants. 
As a result of this definition, we have 



•ik"™ 


( 2 - 12 ) 


/here m is the total mass of the fluid system. We note that because of 
Equation (2-11), dXj^ is independent of species i. 

The rate of change of mass of i^^ species due to all the reactions is 
therefore. 


dm|^ 

dt 


k=l 


dX 


ik 


k 

dt 


i =1,2, 


(2-13) 


or, in terms of nqass density 


6 

i 
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dp. ’*• dX, 5p. 

dt ^ Z-J ik dt St 
k=l 


(2-i3a) 


With this relationship the species continuity Equation (2-8) becomes 


Dc. 
1 

Dt 




r Da. 

L-'skDr 

k=l 


(2-14) 


It will be shown later that the time rate of change of progress vari- 
able dXj^/dt is directly related to the temperat re and pressure variations 
in the flow field. 


It is worth noting that all the above equations are extremely general 
in their validity. That is, they hold true for all flow systems including non- 
equilibrium chemical reactions and irreversible transport effects such as 
viscosity, conductivity and diffusivity. 

2.1.2 The Momentum Equation 

By applying the principle of conservation of momentum to a differen- 
tial element in the flow field, one can derive, in vector for n, 

+ V • (pq?) = + S Vi 

i=l 

where is a dyadic product of the velocity vector Q by itself and T. is 
the body force exerted on the i " species per unit mass. The stress 
tensor a is siometimes written in the form 

I 

5^=-pJ+T (2-16) 


7 
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Here, p is the hydrostatic pressure usually defined as - (<r,j + 0-22 + 0jj)/3; 
5 the imil U nsi*r; and T the viscous part of the stress tensor. 

Putting Equation (2-16) into Equation (2-15), we have the well known 
Navier-Stokes equation for compressible fluid flow, 

• (pqq) = - Vp + V» T + 2^ p.f. (2-17) 

i=l 

In this equation the five terms are, respectively, the non-stationary and 
convective rate of change of momentxim per unit volume, the net hydro- 
static pressure force and the viscous stress force acting on the surface of 
the unit volume, and the body forces per unit voltune. 

By using the continuity equation (2-7), Eqviation (2-17) may be 
written as 

= ^ + 5.75 =.ivp + iv.f+i£p,rj ( 2 - 18 ) 

i=l 

For isotropic, Newtonian fluids, the viscous stress tensor r in Equations 
(2-17) and (2-18) can be shown to be related to the velocity gradients 
(strain rates) by the following formula 

T = (7 • q) I + ^^^2 ^ m)^)/2 ( 2 - 19 ) 

T 

In which (V^) denotes the transpose of the tensor 7q, and and 1/^2 
called the first and second coefficient of viscosity of the fluid, respectively. 
These coefficients are usually temperature dependent. 

When t le body forces are negligible, as is usually the case for 
compressible flow, and if the viscosity effects are small, then we arrive 


8 
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at the familiar Euler equation of motion 

tq-vq — jVP (2-20) 

We note at this point that chennical reactions do not alter the forms 
of global continuity equation or equations of motion. 

2.1.3 The Energy Equation 

The most general form of energy equaMon, according to the first law 
of thermodynamics or the law of conservation of energy, can be written as 

^ Jp(e + iq^)] + 7- [pq(e + ^q^)j=-V .Q + V(a*q) 

+ 2Pi(q + q()-fi (2-21) 

i=l 

where e is the internal energy per unit mass of the gas mixture; q the magni- 
tude of the velocity; and Q the conduction heat flux vector. 


Equation (2-21) simply states that, based on unit volume of a fluid 
element, the total rale of increase of internal and kinetic energy due to local 
and convective changes (the first two term.i) is equal to the heat conducted 


into the volume (the third term) plus the work done on the fluid element due 

I 

to stresses (the fourth term) and the work don^ on the element by the body 
forces (the last term). j* I 


In writing Equation (2-21), all that has lieen neglected is the energy 

l' 

transferred due to radiation, which if necessairy, ' can be added to the equation 
as an exti a tet m. i ; ' i ; I 


N^dless to say, the microscopic quantui|i effects (interchange of 
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energy and mass) are not considered in Ev^uation (2-21) as it has been im- 
plied by the classical first law of thermodynamics. 

Using Situation (2-16) and introducing the enthalpy per unit mass 

h = e p/p 

Equation (2-21) can be easily rearranged to give 

^ P(h+iq^) +V • [pq(h+|q^)j = +7*(T*q-Q) 

+ ^Pi(q + qp • I (2-22) 

i=l 

Or, if we subtract the continuity Equation (2-7) from Equation (2-22), we 
obtain 


^(h+iq^ = i^+}v-(f.q-Q)+j;^Pi(qtq')- (2-23) 




This is the energy equation in which all the transport properties - viscosity, 
conductivity and diffusivity, have been taken into consideration. 


2.2.4 The Entropy Equation 


The gp neral thermodynamics equation with chemical reactions may 


be writte^ 


1 


.( 

i' ! 

J: ; 


Tds = dh - ^ dp - > p;dc. 

p ^ 


' I 


(2-24) 


where s k tl.u entropy per unit mass oftheimixtu^ and the chemical 
potential ^er mole of the chemical si>ecie8' i. jl | ' 


10 
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For a flow system, we can write Equation (2-24) in terms of 
substantial derivatives 


Dc. 


^_Dh yr ^ ZZl 

^ Dt ~ Dt p Dt Z-> i Dt 

i=l 


(2-25) 


Now, multiplying the momentum Equation (2-18) by q, one gets 

n 


q 


1 1^1 ^ 
q . (q «Vq) = - - Vp • q+ ^ (7 • T) . q + p ( 


2-26) 


If this equation is subtracted from Equation (2-23) the result is 


§r-?^=7^:^5-iv.5ti£p.q'.f. (2-27) 


i=l 


where tl» symbol J is a double dot product of two second order terms- It 
is defined as 


3 3 

EZ 

j=l k=l 


- EE Ajk»jk 


80^ 

= E Z ’■jk 

j=l k=l 1 ^*k 

!■ : 

Now, substituting Equation (2-27) into Equation! (2-25), we readily obtain 

I; I 

Tsf ' E>i9i'- VE^ ^ I 

. i=l : i=l j 

I I s 

This equation states that the total s*ts| of entropi production, in energy per 
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unit mass per unit time, is equal to the sum of the individual rate of entropy 
production due to viscosity (the second term}, conductivify(third term), diffusi- 
vity (fourth term) and chemical reactions (last term). The fourth terms on the 
right-hand side of Equation (2-28) are all positive quantities. They are 
sometimes known as the dissipation terms. 


The last term in Equation (2-28) may, with the help of Equation (2-14), 
be written as 


n 


•E 

i=l 




Dc. 

1 

Dt 


--E^ 

i=l 


r DA. , 

E-ik 

k=l 


(2-29) 


which, if the diffusion velocity is neglected, becomes 

n Dc. _r_ DA, 


■ E"i^ = -E^°k m 

i=l k=l 


(2-30) 


n 


with 






(2-31) 


th 


It is to be noted here that the Gibbs free -energy change due to k reaction, 
.AGj^, vanishes for equilibrium chemical reactions.^ 


2.2 THE EQUATION SYSTEMS FOR STEADY NON -EQUILIBRIUM FLOWS 
OF REACTING GAS MIXTURES WITHOUT TRANSPORT EFFECTS 


For steady, adiabatic, inviscid and non-diffusive flows, the general 
equations derived in the last section may be re<^uced to relatively simple forms. 

2.2.1 tL 


e Continuity Equation f: i 

f.i 

Tlie global continuity equation (2-7) becomes 


7*cpq) h oij 


(2-32 
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The species continuity equation (2-14) reduces to 
Dc. r DX 

-dT ' “ 

k=l 


In the following we shall show that the variation of progress variable 
Xj^ with respect to time is directly relat* d to the temperature and pressure 
variations. 

f2l 

Following the method of Kirkwood and Crawford we separate 
into t\'o parts 

Xk = + 4 (2-34) 


The parameter X^ is the equilih-ium value of the progress variable for the 
th ^ 

k reaction for the instantaneous local temperature, pressure and com- 
position. The variable A is the measure of the deviation of the chemical 

tH 

reaction from equilibrium; i.e., it is the lag of the k reaction in its attempt 
to maintain equilibrium. 


The equilibrium value of X^ is determined by the condition that the 
Gibbs free -energy change for reaction k vanishes, i.e.. 


= 0 


or, by Eouation (2-31) 


Since 


at 




■ I 


2.^ *xk Dt "i®- 
k=i : : 


dp j = « SjdT + 


“ iti 


n 




(2-35) 


(2-36) 

i 


(2-37h 
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Equation (2-36) becomes 
n n 


dll ■ Dc^ 

EE dc I Dt 

1=1 i=l ^ 


= As, 


51 .JL\ 

v^A 


5 £ 

Dt 


(2-38) 


where 


^». = E-'iKV 


i=l 


i=l 


are respectively the change of specific entropy and specific volume due to 
uth 

k reaction. 


Now if we substitute the species continuity Equation (2-33)for X. = 

J 3 


into (2-7) we obtain 


i i 1 DT 
‘'ik ‘'fj ^ ~Dt " ‘^^k Di 
i=lf=lj=l ' 




5 e 

Dt 


(2-39) 


k = 1,2 r 

This equation holds true for flow systems of reacting mixtures under com- 
plete chemical equilibrium. 

2.2.2 The Momentum Equation 


Here, we obtain from Euler's equation |Of motion (2-20), 


2.2.3 T 


q . 7q = - 1 I 


(2-40) 


he Energy Equation 


I !l 


Ffom Equation (2-23) the steady state i|e|gy equation for inviscid. 


14 
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adiabatic flows with negligible body forces is simply 

h+ = h = const. (2-41) 

where h^ is the stagnation enthalpy per unit mass of the mixture. The local 
specific enthalpy of the mixture is to be evaluated by 

n 

h s (2- 41a) 

i=l 

in which h^ is the partial specific enthalpy. 

2.2.4 The Entropy Equation 

From Equation (2-28) and (2-30) we have simply 

k=l 


where was given by Equation (2-31). 


In the above, we can cov a total of n + r + 4 equations, namely, n 
equations of (2-33), r equations of (2-39), Equations (2-32), (2 40),(2-41)and 
(2-42). Yet we have a total number of n + r + 6 dependent variables, namely, 

c j , c^, Xg, .... X 2 . p. T, p, h. s, and q. Thus, two more equations are 

needed. The first one has already been implied by the definition of c.. i.e., 
(Equation |(2 -4). 


n 


c. * 1 
1 


(2-43) 


i»l 


The second one is provided by the eqxiation of state of the gas mixture 


P = 



(2-44) 


i 
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where R is the universal gas constant and W. the molecular weight of the 

th U » 1 s 

I species. 

This completes the system of the basic governing equations for any 
reacting flow without transport effects. . 

When the reactions between components of the gas mixture are 
known and the boundary conditions adequately specified, one should be able 
tc solve the non -equilibrium flow system. 
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Section 3 

FROZEN AND EQUILIBRIUM FLOWS 
3.1 FROZEN FLOW 

Under many circumstances, chemical reactions may proceed so 
slowly that there is hardly any change in composition of the fluid during the 
time for il^^ traverse the whole region of flow field of interest. In the 
limit the nnxture may be considered as having fixed composition. Such a 
flow field is usually described as a "frozen flow." Under a frozen flow 
condition, th.; mixture behaves as a single perfect gas. The species con- 
tinuity equation does not come into picture while the equation of state becomes 

p = R pT/W 

u ' m 


where W is the fixed molecular weight of the mixture, 
m 

Under frozen flow assumptions, the basic governing equations for 
steady flow of gas mixture without transport effects are as follows: 

Continuity: 



mentum: 


7 • (P q) s 

i i 


q . - i 

p I . ' 





(3-1) 


(3-2) 


(3-3) 


Ehcrgy: 



L.MSC-HREC Dl62220-m 


Here, h is a function of temperature only. 

Entropy; 

s = constant along a streamline (3-4) 

State; 

p = 

The above equations are the usual gas dynamic equ i for 
inviscid, non-conducting, non-diffusive compressible flow without 
chemical reactions. 

It should be stated here that the validity of Equation (3-4) rests upon 
the assumption of zero transport effects. This means that ciiscontinuities 
such as shock waves should, strictly speaking, be excluded, since the 
appearance of shock waves indicates the existence of transport phenomena. 
However, experience tells us that shock waves usually occur in a very 
narrow region (in this region viscosity and diffusivity predominate) and we 
can still apply the isentropic criteria on each side of the shock wave but not 
across it. 


EQUILIBRIUM FLOW 


Equilibrium flows are defined as flows in which the mixture is at 
chemical equilibrium at all instants and ever^here in the flow field. Phy- 
sically, this is approximately the case when t|ie chemical reactions in the 
mixture proceed so fast such that the reactions are "completed" in a time 
interval much smaller than the time interval we hav* 'eeson to be inte- 
rested in, A reaction is said to be "completed^ when there is no further 
apprecia >le change in the composition of the mipeture. 


8 
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In this case, the set of governing equaf’''ns may be again sin^-Uified 

(•o give: 


Global Continuity: 


V.(pq) = 0 (3-6) 

Momentum: 

q • Vq = ' ^ VP (3-7) 

Energy: 

+ I (3-8) 

where h^ is the local equilibrium value of specific enthalpy- 
Entropy: 


s s constant along streamlines 


(3-9) 


This is a consequence of Equation (2-28) and the remarks following Equation 
(2-31). We note that the flow is again isentropic in ® shock-free field just 
as in the case of frozen flow. 


We note at this point that chemical reactions do not =lter the forms 
of the above equations as compared to frozen f)ow case. 


Swte' 




(3-10) 


19 
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Section 4 

UNCOUPLING THE FLOW AND THERMODYNAMICS FOR 
STEADY. FROZEN OR EQUILIBRIUM FI OWS 


For steady flow equation (3-13) can be written 


n n r 


EES "ik ■>j 5^7 " 

i=li=lj=l * ^ 


0 (4-1) 


Since all terms except the gradients of Ay T.p are scalar quantities we may 
write 


n n r 9M • „ 

^• EEE "ik'ij 

i=li=ij=l < 


k ^(i) JP = 0 


Thus the thermochemical behavior does not depend on the magnitude of the 
velocity, and assuming that it does not depend on the direction, we have 


Ik •'ij ac. 
i=ii=ij=i ^ 


Wx°-<>s.VT+ A i) Vp = 0 


(4-2) 


Examining now the behavior of the thermodynamic system over an infini- 
tisimal step dL 


n n r 


EEE ■'ik -ij 5^7 - ^'k-*T + 

i=li=lj=l ' ' , 


A(-)^dp = 0 


(4-3) 


Obvioush’ the above equation does not depend on position or velocity and 
hent.e is entirely uncoupled from the flow problem. Using the same approach 
on Equati >n (3-11) we find; , : ; 


dc. - V.. dX? 4 ;() 

i I ik k , , 

k=l , fl , 


(4 4) 
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The thermochemistry equations (3«11) and (3-13) have been completely 
divorced from the flow problem and describe the behavior of a mixture of 
perfect gases under infinitesimal changes in pressure and temperature. 

It is no longer necessary to use Equations (4-3) and (4-4) in their 

differential form. Systems of equations describing the behavior of reacting 

mixtures are well known and many sophisticated computer programs exist 

to facilitate the solution of these equations. The NASA-Lew»is thermochemi- 
[3] 

cal progranr •'was chosen for use in this study for a variety of reasons and, 
due to the fact that special features of this program were utilized in simpli- 
fying the subsequent flow field analysis, other thermochemical programs 
may not be immediately compatible. This need not, however, preclude 
their use. 
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Section 5 

PROPERTIES OF THE WORKING GAS 

Given the initial reactants, combustion pressure, and total energy 
the thermochemistry program chosen for use in this analysis evaluates, 
among other things, the following quantities; 

T. s, w ,r 
m 


where >', the isentropic exponent, is defined by 



(5-1) 


At constant entropy, a series of lower pressures are then evaluated. Aieach 
step the enthalpy of the mixture is computed in addition to the quantities 
listed above* From the conservation of energy, the thermochemical analy- 
sis is related to the flow systems by calculating the energy change from the 
original combustion value* In the flow problem the energy difference is the 
kinetic energy of the gases. A velocity may then be fo xnd corresponding to 
each pressure. In order to describe the effects of shock waves and other 


irreversible phenomena the cnamber pressure is decreased and the process 
repealed. In this fashion a tabular description of the behavior of the thermo- 
chemical system is constructed. Entropy and^velocity are chosen to be the 
independ|nt quantities and specification of the^e parameters is then suffi- 
cient to Uniquely describe the gas properties, f: | 


Computer economics dictate that the tabular description occupy as 
little spa :e as possible. Because of this an accurate interpolation scheme 
is necessary. Before attempting an explariation of the interpolation scheme 
a definition of local reference conditions is in.c^der. Usually, for an ideal 

» , i .1 ;i : 
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j*as, reference conditions are chosen to be iscntrupic stagnation conditions 
To arrive at local conditions one need only employ well known gas dyr rr.'c 
equations. For a gas whose isentropic exponent and molecular weight is 
varying, however, the situation is more complicated. A reference con- 
dition is defined which is produced by iseatropically stagnating an ideal gas 
whose isentropic exponent and molecular weight are equal to the local mix- 
ture value. The reference temperature is 

T„-T.|^q^ (5-2) 


and the reference pressure is 


r 



where R = R /W 

u m 


(5-3) 


To interpolate between two velocity values (1, 2) use 

R = (1-h^.) R^ + h^.R^ (5-4) 


where 


y = >'2 



(5-5) 


(5-6) 


(5-7) 
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and to interpolate between two entropy values use 


where 


R = (1-hs) Rj + hsR 2 


(5-8) 


y = fi-hs) Xj + hs 


(5-9) 


p = Pj exp 


/ . ^2 
hs xn-s— 


( 5 - 10 ) 


T = (l-hs) Tj + hs 


( 5 - 11 ) 


bs = 


be s 


S-8, 


® 2*®1 


s/R-8j/Rj 

®2^^2“®l'^^l 


In the foaowing development of the flow et^uations no further mention is 


made of 


he nature of the thermochemistry of the 


of the prpperties is assumed to exis^ so t^t it 


gas is ideal, froaen, or in chemical equilibriilxhi 


U 


is immaterial whether the 


gas* A tabular descriptionj 


U ' 
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Section 6 

THE METHOD OF CHARACTERISTICS SOLUTION 
6.1 DEVELOPMENT OF THE CHARACTERISTIC EQUATIONS 

The continuity equation or conservation of mass is, in vector form, 

V MPq) = 0 (6-1) 

while the conservation of momentum is 

- i Vp = Vq^2 - qx (Vxq) (6-2) 

The thermodynamic relation is 

T ds = dh - ^ 

P 
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Expansion of (6*1) yields 


7 • (pq) = P7 • q + q • Vp = o 


(t-6) 


which when combined with Equation (6-4) gives 


pV»q + q* 


,oP 


a 




But since q »Vs = 0 aud is defined to be the speed of sound we have 


\7 • q “ ® 


(6-7) 


Recalling the definitio* equation (5-1) 






so that 


a*" s rRT 


(6-8) 


Choosing now a cylindrical coordinate system (Figure l),in which the flow 
angle $ is defined to be the angle between the^ velocity vector and the x axis 

and where r is the distance from that axis, and in this system expanding 

‘ l4i 

Equation (6-7), we have for axisymmetric or ,|twO dimensional flow.^J 


('■?)**•(■■ 


2 

V 

“I 

a 




9r 


2 [ax dt 


(6-9) 


where <rks 0, 1 for two dimensional ^r axis,. ^ 
aider a fine in the r.x plane inclinet^ to thp x ifi 
the express ion for this line (normal 

■ ^ 

3^ 



efric flow respectively* Con 
by the angle tan"^^, 
hLaracteristic line) is 


II 


then 
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But along this line 


and 


du = |!id:< tfi dr 
ax dr 



dv = 




dx 


then 


au du „ au 
dx ■ dx " ^ ar 


and 



I 


1 


/ 



( 6 - 10 ) 


( 6 - 11 ) 


|8 
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Substitution of Equations (6>10) and (6-11) into Equation (6-?) yields 
simplification, 



Recalling Equation (6-5) and rearranging yields 

2 

7s s qx(7xq) 


where, for two dimensional or axi symmetric flow, 

jdx arj 


and 


a^ ds 
m dr 


= u 


gv gu 
9x ' gr 


but 


where n 


gs _ gs gx gs 
dn ” gx dn gr! gn 


is normal to the streamline and 


dx 
gn I 








after 


(6-12) 
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and 


dji 

dn 


u 

q 


so that 


2 

go ■ a ds 
dx ~ dr yRq do. 


Substituting of the above result into Equation (6>12) yields 


I, XX \ du . 1 L v^\ dv , o*v 

"■7 3 ^ *-71 3;*- 


'IV . 1 I, V 

r-Ti 


as 

yRq dn 


(6-13) 


2uv 




1 

0 


a 


^ = 0 
dr 


It is desirable to determine whether there exists some /3 for which the 
coefficient of ^ is zero. 
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Lot 


u s q coB$ 

V s q sin0 
ft u sin“^(a/q) 

Then, after much manipulation, 


s tan ( 0 + ft ) 


(6-14) 


Combination of (6-13) with (6-14) yields then 


2l 


UV 


+ cot(0 + ft)ll- — ^ 
a V a^ 


2 ^ 
_a ^ 

XRq dn 


Now 


ds = dt = ||dn = |i sinftdf 


where t is taken along the streamline and d| is measured along the character- 
istic line. But 


dl = 1 + 


m 


1 


dx s sec(d7ft)dx 


Then finally after much .algebra 

= 0 (6-I5) 


Equation|(6-15) is known as the com] 
variatioijin flow properties along 

dr 


tai^ (^f|fi 


ipatabjlity ;e|qi|atlon and describes the 

(6-16) 
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These lines are the so-called physical characteristic lines. Locally they 
are tangent to the Mach lines. Equations (6-1 S) are often called the hodo- 
"raph characteristics for two dimensional flow. In "cisymmetric flow the 
existence of the physical dimension causes a coupling of the characteristics 
so that Equation (6-15) cannot truly be referred to as holograph character- 
istics. 

It is interesting to note that the only reference to gas properties in the 
derivation was the use of the definition given in Equation (5-1). 

Examination of the four equations (6-15) and (6-16) reveals that there 
are five unknowns. An additional relation is provided by as summing a linear 
variation in entropy between the known data points. 

6.2 FLUTE DIFFERENCE SOLUTION OF THE CHARACTERISTIC 
EQUATIONS 

In order to solve the differ^ itial equations (6-15) and (6-16) it is first 
necessary to v'rite the equations in finite difference form. At times in the 
flow field certain conditions are known which allow some of the equations to 
be discarded. These are, of course, the boundary conditions. 
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■ ^1,2 

^1,2 ^1,2 

and Stakes on the values 

I 

1 

1 

for interior solution 

s, = < 

• 

for lower boundary solution 

1 

* 0 

i 

for upper boundary solution 

while 

1 

1 1 

for interior solution 

^2 “ 

) 0 

fci lower boundary soli ' » 


1 1 

for upper boundary solution 


which correspond to the Figures (2a) - (2c) given below 
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Figure (3c) 


Notice f’hat for the boundary condition solutions the Equations (6-17) be- 
come the streamline equation when S is zeroJ For a solid wall solution 
the flow jangle is known while for a pressure [boundary the pressure is 
known. 


Five basic cases will be discussed. Thede ate: 


a. interior point solution} 

b. upper wall I I 

c. upper free boiuv 

d. lower wall 

e. lower free bourn 
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The compatability equations are: 


a. interior point 
*^3 = 


" ^2 * Qj qj + Q2^2 ^1 ^2 ‘ ’ ®2 

Q1+Q2 


^3 = ^2 + Q2<‘l3-V-°2 + ®2 


b. upper wall 


- ^2 ^2 " ®2 ®2 *^2 
^3 5 : 


(6-19) 


(6-20) 


(6-21) 


where 0^ is given by the wall equation 
c. upper free boundary 


«3 = «j + Q^(,3-,3).G2+B2 


and q^ is found from the local pressure 


d. lower wall solution 


* 


- ^3 + Gi - + Qj qj 


where $l is given by the wall equation 
lower pressure boundary 


^ a 

I ! 


and q^ il given by the local pressure 


35 


■ ill 


i I ■' 


( 6 - 22 ) 


(6-23) 


1 ' 

> ( 
1 

- Bl 

(6-24) 

i 

1 

1 

! i 
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where 


cot/! 


1.2 


‘ 1.2 


' 1.2 


1,2 


**1.2 

o-sia^^ 2 
'1.2 
sin/! ^ 2 M 


1.2 


1 ,2 ^^3 ~ ^ 1 . 2 ^ 


^1.2 ^1,2 


(6-25) 


and is given by 


a. interior point 


So = 


4. fillLlfil 

F. + F, 


1 


(6-26) 


b. and c. upper boundary 


‘3 = «1 


(6-27) 


d* and e* lower boundary 


83 = 82 


(6-28) 


while 


point. 


1.2 


siny , 3 (xj - X| 


IfZi 

f I 


( 6 - 29 ) 


An iterative solution is employed to determine the properties at the new 


'or the first pass of this solution the barred values are approxi* 


mated m the conditions at the two known points. 


I 

. t 

d. 1= e 


I 


1 .1 
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for the first iteration. After the aporopriate set of equations has been 
solved a new estimate of the barred v^alues is made. For example, 


h 


2 


The iterative process is continued until the desired convergence is reached. 
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Section 7 

THE OBLIQUE SHOCK SOLUTION 
7.1 DEVELOPMENT OF OBLIQUE SHOCK RELATIONS 

Figure 3 illustrates a stream tube passing through an oblique shock 
wave. This wave, which is extremely thin, will cause an almost instantan- 
eous rise in pressure and temperature. For some distance downstream of 
the shock wave (in a reacting gas) a non-equilibrium zone will exist followed 
by a return to chemical equilibrixun. The following analysis discusses the 
fluid flow properties in such a way that the non-equilibrium process need 
not be specified in order to arrive at an exact solution for the gas properties. 

Consider a control surface as shown in the figure. The conservation 
of mass yields: 


■ ^1 *^ 1^1 ^ ® 


(7-1) 


Conservation of momentum gives 

Pi ■^1 

MA k * * ST 


— / Pz^2^ 


"(Pl + Pj qj) Aj fn + (p^ + Pj W - 


f ■ /^Pn.E. 

^.E. ^N.E. 


dA = 0 


Since each 


streamline locaP / undergoes the same process the last two terms 


of the abov4 equation are equal and opposite. Therefore . 


a . 


A, in' 


i 1 ;l : '■ i 
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So that, after substitution of (7-3) and (7-4) and setting each component to 
zero. Equation (7-2) becomes, 


2 2 
p2q2Cos(€ - S)A2 - p^q^A^cos^ = 0 

(p, t P, qf) A, s.nt t cos( -fpj + P^qfjA^smie-S) 


(7-5) 


p2A2Cos(€-8) 

tan(£-$) 


= 0 


(7-6) 


But from geometry it can be seen that 


2 _ sin(g-S) 


sine 


(7-7) 


After substution of Equation (7-7); (7-1), (7-5) and (7-6) become 


^2^2 - Pi 'll siu€ = 0 


(7-8) 


P 2 <l2 ' Pi = 0 


(7-9) 


?2 + Pz^2 - Pi - Pi sin^« = 0 (7-10) 


The above set of relations contains €, 5. p »P 2 * ^2 unknown quantities, but 


?2 = PC^z-^lz^ i Pz = P<»2''*2> 


(7-11) 


So that if on<i variable, say e. is taken as an independent parameter the 
remaining unknowns i$, q 2 ,S 2 )>n^y be found by an iterative solution. These 
equations ar s, of course, formally the ^ame as the ideal gas solution. The 
difference lins only in the variation of pressure etc., with entropy and velocity. 

I I li i -i I 
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It is impossible to determine the location of the new equilibrium shock 
point location without a detailed description of the non -equilibrium reaction 
process. It will be assumed therefore that this zone is thin and that no signifi- 
cant errors are introduced by letting the downstream physical location lie 
on the upstream location. 

It has been pointed out previously that the characteristic equation 
derivation was based on neglecting transport properties and as such is 
necessarily restricted to continuous regions only. The oblique shock wave 
relations derived here then are patching lines betw’een the continuous regions. 

7.2 ITERATIVE SOLUTION OF THE OBLIQUE SHOCK RELATIONS 




41 
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From calculus 


ao. dC. 


Now 


dCy dG- 

‘*°2 ' « 5 T '*'>2 557 ^‘z 


(7-15) 


5G 


1 


"l .2. 


^2 = 


("“' 2 ) 


= af; ('”'’ 2 ) 


^®2 H‘^2/ VzI 


^ = P2 8^(^''P2) -fe) ‘'f “‘'■^‘P2^('"P2) 


M 2 


9G 


557 = P2 efr ('"P2) - (k) P2 6^ <'’'P2> 


Rather than calculate the partial derivatives munerically by perturbing the 
iunctions inp 2 > approximate values for these derivaties will be found 

by assuming that locally the gas behaves ideally, that is to say 


dRr 


ds^ ” ^ 


^R, dTo, axo- 


ds. 


aq, 


8PO2 


42 
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ds~ ~ dq. 


:: 0 
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so that 


0{lnp^) ^ d(tnp^) ^ j 

OSo ^ 


(7-17) 


and 


Pz> S y _2_ 


dq. 


^2 dq2 ~ ' 


P2 


(7-18) 


writing (7-15) in finite difference form: 


a/-> 1*^) « 

o <“> = ^ L 

1 I aq2 


dG 


(n) 


(n-^n.o 


ds-» \ 2 




(n+l) (n)] 

2 ‘®2 


Since the root Gj = G 2 = 0 is desired, G 2 ^*^^^^ are set to zero, 

resulting in 


fiC' dG 

(nt ) ^ (n) ( (B) . o (B) I 

2 aq2 1 





dG2^”^ dG.^"^^ 


as 


2 , ^2 

i 


as, aq 


inM 

t) 


(7-19) 


43 



LMSC-HREC DI62220-IU 


and 


q(n+l) , ^(n) 


as. 


("1 


.(n+1) fn) 


f6G 
1 

c>q, 


(n) 


(7-20) 


The iterative solution using Ecuations (7-19) <and (7-20) is continued until 
the desired convergence of and G^, is reached. The solution is com- 
pleted by 


*. ■ "1 
S ~ € - sin 




sins 


(7-21) 


The first guess to start the solution is an ideal gas solution to the set of 
equations. If it is indeed an ideal gas under analysis the first guess is exact. 
These relations are 


S ~ ~ tan 


-I 


tans 


1 




,, 2 2^ 
Mj Sins 


Xj + l 


= 


coss 

^1 cos(S-5) 


(7-22) 


^1 

in 

2 XjMj^ sin^S - (Xj - 1 ) 

+ Xjin 

tan(s-S) 


®2 ■ ®i 

Xj + i 

tan€ 






* 
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Section 8 

EXPANSION CORNER - PRANDTL -MEYER FAN 

In some cases the flo'.v may be required to negotiate a sharp expansion 
turn. The problem becomes two dimensional at a sharp corner (it is im- 
possible to conceive of an expansion corner on an axis of symmetry) and 
may be treated with a Prandtl-Meyer expansion. 


q + dq 



Figure 4 


Since a Mach wave will support pressure changes only in a direction 
normal to itself/”^-^ 


dv = qd^ 
du = dq 


45 
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or 


du 

• 5 — = tanu 

dr 


1 

^ M ^-1 


dd 


= ^M^-l 


q 


( 8 - 1 ) 


The solution to Equation (8-1) is a straightforward numerical in^e 
gration for the case of a known final velocity (free -boundar y case). If the 
turning angle is known, however, and the final velocity is not known, an 
iterative solution is necessary to determine the upper limit. 



(8-2) 


In the mesh construction to be discussed later a fan of rays must be 
generated to allow a numerical description through a large turning angle. 
The turning angle is subdivided into a number of small turns, each of which 
is in^ grated numerically. Corresponding to each o. 'mr.ll turns is a 

Mach wave or characteristic line emanating from the co. 
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Section 9 

SHOCK NEIGHBORHOOD SOLUTION 


Due to the fact that the shocic calculation will be employed under 
several different flow conditions a general setup and notation is used. 
There are six basic types of calculations as shown in Figures (5a - 5f ). 



Figure 5a - Interior Right-Running Shock Wave 



Figure 5b - Interior Left-Running Shock Wave 


47 
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Figure 5c - Solid or Free Upper Wall Interac. '.g 
with Left-Running Shock 



Figure 5d - Solid or Free Lower Wall Interacting 

{ with Incident Right-Running ShocK Wave 
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Figure 5e - Solid or Free Upper Wall Interacting with 
Reflected Shock Wave or Attached Shock 
Wave with Insufficient Downstream 
Information 



Figure 5f - {Solid or Free Lower Wall Interacting with 
Reflected Shock Wave or Attached Shock 
jWave with Insufficient Dow'nstream 
Information 
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Let 

8 - = 1 . - I , - 1 . 1 . 1,-1 
8 ^ = 1 . 1 . 0 . 0 . 1. 1 
6w= 1.1. 1.1. 0*0 

for the six cases, respectively. 

The angle that line 6-8 makes with the axis is 

while the angle that line 7, 9 makes with the axis is 
while the angle of the shock wave is 



(9-1) 


For an initial approximation let ^8^ + ~ ^2 ^ *2 Compute the physi- 
cal location of which is just the intersection of the shock wave and 

the line 6-8. 

A linear interpolation between the flow values at points 6 and 8 is used 


to determine the local flow properties at point 4. A shock angle is then 
determined. i i 
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The oblique sHocic solution of section c is used to determine the flow 
properties at 3 across the shock wave. 

It is now necessary to make a characteristic cal ulation to determine 
whether the new shock solution is the correct one. Li order to do this a 
physical characteristic line is drawn from point 3 and the intersection of 
this line and line 7-9 is point 5. In the cases shown in Figure 5e and 5f 
point 9 is not known. An estimate compatible with the shock slope is made 
using a characteristic solution. 


For a first approximation to the slope of the line let 

(n.o) ^ 0 (n) ^ j 

5 3 and let fi' ’ ' = n' ' 


then the characteristic line makes the angle 


(n.m) _ 




(9-2) 


The doublet may now be found and an interpolation between the pro- 

perties at points 7 and 9 yield the estimate of the flow properties at 5. This 
information is used to improve the estimate of in Equation (9-2). The 
process is continued until convergence is reached. The physical character- 
istic line passing through point 5 terminates at the shock at point 3. If the 
compatability equation is now solved along this line a measure of the in- 
accuracy of the shock slope will result. Using all the known information 
concerning points 5 and 3 solve for the flow angle at 3. 

= V"' - Q (^3*"’ - G - Si B (9-3) 


51 



LMSC-HREC D1 62220 -m 


where 


Q = ^ 


^ sin^ F 


sinU cos fi (s^-Sg) 
R ? 


sinfi (x*-xj 
F = :i - - - 


aad, of course, the barred values are averages over the step length. 


An error may then be computed; 


but 



and the solution is perturbed until the root is bracketed. The solution is 
driven to convergence by the method of false position. 

It should be pointed out here that the use of linear interpolation is a 

i 

matter of convenience only. It would be possible, but in general far from 
practical, to generate by exact means those {joints which have been approxi- 
mated. 
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Section 10 

MESH CONSTRUCTION FOR INTERNAL FLOW 

The calculations described previously are point or small region 
solutions. Some process must be de/* cd which successively employs 
the proper calculation at the proper time in order to describe the entire 
field* In order to facilitate a description of the mesh construction process 
let represent the total knowledge of flow properties at a point in the field. 
Also let the expression 


♦ = V 


stand for properties at a new point which are computed as a function ^ of 
(m) other points. There will be basically six such functional operations 
^ stand for input point, interior point, bound- 

ary point, attached shock point, shock, and Prandt-Meyer pcin«.s. In addi- 
tion the superscript (u) will indicate that the operation is to be performed 
in the presence of an upper boundary while (L) indicates a lower boundary. 


Due to the complexity of handling multiple shock waves, a single 
shock wave restriction will be imposed. This shock wave is arbitrarily 
chosen to be of the right running family. This type of problem will :ost 
frequently occur in cases of interest. If a left running shock wave occurs 
the problem is simply inverted. 

i 

The choice of ri^ht running shock waves also dictates that left running 
characteristic lines be followed in the calculation. This allows one to re- 
tain a minimum of information, i. e., a known characteristic line (hereafter 

referred to as (j)) and a line in the process qf Ibeing computed (k). 

I ) I 
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To begin the problem all necessary boundary conditions must, of 
course, be supplied. In addition a starting line containing N points which 
are designated (n = 1, . . . N) must be supplied. 

Figure 6 illustrates a flow field in which there are no discontinuities 
and in which the mesh construction is terminated when the region of interest 
has been computed. 



In region I the left running characteristic lines initiate as input points 
and the mesh construction may be described by; 







i = 1 

i = 2 ... (2n-2) (10-1) 

i - 2n-l 
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where n varies from 1 to N and <{>. . represents the flow properties at the 
th ^ 

i point on the k line. For instance, in calculating the fourth point on the 
fourth line shown in the figure, line three is known in its entirety and line 
four up to and including the third point is known. The above set of relations 
says that point three on line four and point three on line three will deter- 
mine, through the interior point solution, the next point (four) on line four. 


For region II we have; 



“ 2....(2N-2) (10-2) 




As a new line becomes completely defined it may be referred to as j 
and the process continued indefinitely. 


It is possible to combine regions (1) and (2) into a more general 
scheme if a variable is defined which takes on the value (1) in region I 
and (0) in region n. At this time the number of points on the j line (i,p ) 

and the number of points on the new line (i_ ) are defined. Then; 

*k 



but 


^*‘^ik,k''^ij. ij = 


'^^^ik.k'^ij, j* iJ * 



i.y 


i-1 




1 

2,...i,_ -1 (10-3) 

^k 
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Obviously would have been initialized to (-1) prior to the start of the 

j 

calculation. When the line is finished i„ is set to the current value of i_ . 

J k 

Thus the process for computing the entire flow field fc .* such a simpli- 
fied case is described by the set of expressions (10-3). In general, however, 
discontinuities will arise so that a more flexible description is necessary. 

If, by some process, points were discarded from the (j) and (k) arrays and 

the number of points lost is i^ and ig respectively then (10-3) becomes; 

j k 



*2.j> 

J » 

ip'* /A A . ik = i-1 ^ — i 

ik, k' ij. j' ij = i-2ij^ + ig +ig^ > 


1 

2 ... 



i -1 (10-4) 


but 






where the tilda over and i^^^ indicates that current values are to be used. 
"C.tese variables are reset to zero at the beginning of each new line. 


To illustrate this, imagine that points (1) and (2) have been computed 

«>o 11 and that after point(3^ had been computed in the normal fashion 

it was necessary to discard it. The next point to be computed would then 

be (3^) but if for somei reason it was necessary to discard point 5 on the 

(j) line then the point 0^) would not exist. Therefore a point has been 

deleted on each line (i^ ~ >1) and the diagram and the set of equations 

®k 

(4) indicate that point (6) on line (j) and point (^) on line (k) would be used 

I 

in the computation of (3) on line (k). Alsq the 'total number of points would 

I : it i 
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have decreased correspondingly. Note that discarding points is the same 
thing as discarding right running characteristic lines. 


It is now possible to include a shock wave into the logic scheme. 
Since it is a mathematical requirement that characteristic lines of the same 
family as the shock are continuously intercepted by it, the ability to discard 
points was necessary. If this was the only mechanism for discarding points 
then the logic process would be; 


1, 

x.k = 


ik=i-l ) 
ij 

‘ =‘»k-‘Sk^* 


^1 ik, k IJ, y 




i= I 


x*2,... X 


(10-5) 


'^I<Vk*‘^ij.j> 


= li = i -i +2 

ij = i-2ij^+l+ig + ig J Sj^ 

^ J 


... .) 


B-'ik.k- ig + ij 

^ j 


XSX„ 


k 


t • • 1 


where 


^k N \ 


and where i is defined in much the same fashion as i„ , which is; 

•k i ^k 


*. ■** ■ * 
k j "i. 

I -1 


where i is the location of the upstream shocjk point on the (j) line. 

•j i Mil 
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Figure 7 illustrates the mesh construction when a shock wave 
is present. 



Figure 7 


In this example ® and i = 8 so that i =7. Also i_ s 16 so 
W ®j ®k 

that i— would normally also be 16. The (k) line is computed up to point 7 

and the shock solution is then employed. In this case the shock solution 

finds that three points of the (k) line fall downstream of the shock (minimum 

is one) while two right running lines (points 10 and 11 on the (j) line) also 

are intercepted by the shock wave. Thus in this example i^ s 2 and i» s 2. 

k j 

The set of equations (10>5) then says that the double shock point should be 

points 5 and 6 on the (k) line and that the total number of points on the (k) 

line has decreased to|12. Note also that the value of i to begin the next 

i ®i 

line must be change to 5. '' 


$o far no mention of how the shock wav^ begins has been made. 

^ r j 

There are two types of shock waves considejred; the attached shock wave 
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which arises due to the flow being forced to negotiate a compression corner 
on the upper boundary, and the envelope shock. The first of these is easily 
detected from the boundary conditions and is initially of finite strength. The 
second type .is detected by a mathematical discontinuity in the mesh con- 
struction (crossing of right running lines) and is initially of zero strength 
i.e., a Mach wave. An example of the compression corner solution is given 
in Figure 8. 



Figure 8 


The computation of the (k) line is completed without any prior know- 
ledge that a compression corner exists. A check is made after the bound- 
ary solution and the boundary information indicates that a compression 
corner must be treated. A linear interpolation is performed between the 
boundary point on line, (j) and the fictitious boundary point on line (k) in 
order to determine thf^ flow properties at point u. An oblique shock calcu- 
lation is made where the turning angle is known. Using this point and point 
6 (i.p -1) a new virtual point {i) is computed, i is set to i_ and the 

j I ®k ^k 

shock solution illustrated in Figure ($e) of Section 6 is employed. This 
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shock solution completes the (k) line in the proper fashion. In this example 
i.p = 6 and the next line is computed as previously disc ssed. 

The envelope shock is detected by a crossing of right running char- 
acteristic lines as shown in the figure below. 



Figure 9 


In this example point (5) on the (k) line is found to fall in a previously 
described region (the region between points (3) and (4) on the (k) line). This 
discontinuity in the solution is interpreted as a shock wave. If the grid size 
were chosen small enough the shock wave would initially be of zero strength. 
Point (5) on line (j) is chosen to be a point which lies on the shock wave and 
the shock solution is c.rrplo/ed. The results of this solution are stored in 
the normal fashion on the (k) line. Obviously the only difference between 
this situation ar.ci t^'eiLtment of a previously developed shock wave is to 
modify the (j) line sudh that it appears to the logic scheme as though a 
shock wave crossed the (j) line at point (5). | , 

; I' 

When an expansion corner is encountered on the upper boundary the 

I ; it I 
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logic schcnic is modified locally. Figure 10 -ivuUratcs the mesh con- 
struction in the vicinity of such a corner, 



In this case point (5) on line (k) is expected to be a boundary point. 

It is discovered however, that an expan on corner must be negotiated by 
the (k) line. A point (6) on the (j) line is found by interpolation. A Prandtl- 
Meyer calculation is employed and the fan of points is stored in the (j) line 
above (6). The total number of points to be expected on the (k) line is in- 
creased accordingly and the normal logic scheme will now complete the 
new line. The next line is calculated in the standard fashion. 

An expansion corner on a lower wall is ^somewhat a more complicated 
situation. Since the calculation no longer utilizes the input line, the lower 
wall expansion fan miy be stored in this area. The set of relations, Equa- 
tion (10-5) is modified to that of Equation (10-6). 


I 
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«'i<Vk- *130* ij=r+(v2)v+i 


I = Z, . . . i 


■k 


% <*k- *3’ 


4». 


ij = 


1=1 -lj5 

®k ®k 

l = i -i- + 1 

®k \ 


ik = i-1 


) 


kj) kk k 


'^B^'^ik.k' “^ij, ij = i+{y2)ij^+ + | " 


= Ir 


where 


% ^ ‘t *<^-V*N-‘Sj-‘s^ 


and where 


i 

s 


k 


3 


for which = 1 until all fan points are used up. 
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Section 1 i 

NON-ISOENERGETIC FLOW TREATMENT 


The bulk of this document is concerned with the discussion of the analysis 
of isoenergetic supersonic flow fields. A computer program was written to 
perform the necessary calculations. Subsequently the computer program was 
modified to permit the treatment of non-isoenergetic flowfields. Although this 
treatment is straightforward, decisions were made in che form ox .11 equations 
and coding of the original program which were not the most advantageous 
approach when the non-isoenergetic flow situation was considered. In particular, 
the compatibility equation (6-15) would have been written in the pressure rather 
than the velocity form. In considering the non -isoenergetic analysis the most 
expedient way of modifying the computer program was chosen. The development 
leading to the modified program (consistent with the constraint of minimal impact 
on the coding rather than straightforwardness of analysis) is presented in this 
section. 

To begin, the development of the species continuity (starting with equation 
2-1} could be replaced by atomic conservation equations. Moreover if the 
conservation of those atoms associated with the fuel and the oxidizer were 
considered then, for steady state, 

V • (pj?) = 0 (11-1) 

and 

7 • (p^, q) = 0 (11 -2) 


would result. 
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If the weight flow ratio of oxidizer to fuel (O/F ratio) is denoted by i| 
then (11-1) and (11-2) are satisfied if, 

q. VI? = 0 (11-^) 

and 

V-(pq) = 0 (11-4) 

The assumptions inherent in arriving at equation (2-40) did not involve 
isoenergetic flow so that the momentum equation remains valid. The energy 
eqtiation (2-41), however, remains valid only along a streamline and must be 
replaced by 


h + 1/2 q^ = h^(q) (11-5) 

The equations of Section 3 are modified in an obvious fashion while 
Section 4 remains unchanged. In Section 5, however, a new variable is intro- 
duced into the thermochemistry determination. Since the discussion of Section 
5 is pertinent to a single O/F ratio, the gas properties description must be 
expanded to include a variable O/P. To do this the discussion of Section 5 may 
be followed for each of two values of the O/P ratio which bracket the desired 
O/P and a linear inte tpolation used to generate the required information. 

In Section 6 equations (6-1) and (6-2) are valid but now 


P = P(p.s,i|) so that 


VP = 


' ' P.»? i ' ''P.s 


( 11 - 6 ) 


Expanding (11 -^4) and combining with (11-6) stUi yields (6-7) since q • V* ~ 
q » = 0. In equation (6-7), however, it is understood that a^ = (BP/dp) . 

i., ill 
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The definition of Y culminating in equation (6-8) is similarly qualified. 
Since (6-9) is generated by expansion of (6-7) it remains valid as do all the 
steps culminating in equation (6-12). 

To proceed beyond this point a deviation from the previous approach is 
taken. Expanding the momentum equation (6*2) results in 


1 8P . ^ (Bv 3u\ 

-5?) 


(l!-7) 


p 8r 




( 11 - 8 ) 


but from the chain rule 

o _ 9 8x ■ B dr 

9n ” 9x 9n 9r 9n 

9 _ _^8x 3 9r 

9t ” 9x 9t 9r 9t 

where (n, t) are respectively normal and tangential to the streamline. Hence, 

8x _ v^ 8r^ _ u 

8n ” q ’ 8n q 

It then follo'xrs that 


and that 



l 


(11-9) 
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Notice that the last term in (11-9) replaces 


_a_ 

yRq 3n 


in the previous derivation. Follov'ing the previous derivation with this in mind 
the counterpart of (6-15) is easily shovvn to be 


— n XRT — r ' ' 


The above result is the compatibility equation used in the isoenergetic flow analyst 
Notice that by continued manipulation the pressure form of the compatibility 
equation will result. That this form of the equation is unaltered by the non** 
isoenergetic analysis is not surprising since it is constructed entirely based 
on the momentum equation and the global continuity equation (neither of which 
are altered by the non -isoenergetic flow phenomena). 

The numerical solution to the governing eqtiations is not greatly affected 
by the modification to the compatibility equation. In equation (6- 1 8) 


As 


1,2 



' Pl,2 


" *'>l,2^‘'l.2) 


while in (6, 25) 


( 11 - 11 ) 



a ) 4= _JL_ /^ii.,2 

1.2 Ti2'Pi,2 



( 11 - 12 ) 


Now equations (6-26), (6-27), (6-28) are the finite difference analog to the 
streamline equation. In a similar fashion we may write 



